In [1]:
import json
import numpy as np
import pandas as pd

from keras import regularizers
from keras.models import Sequential
from keras.layers import LSTM
from keras.layers import Dense
from keras.layers import Activation
from keras.layers import Dropout
from keras.layers import Flatten
from keras.optimizers import Adam
import tensorflow as tf
from datetime import datetime 
import shap
import matplotlib
from matplotlib import cm
import plotly.graph_objects as go
import plotly.express as px
from keras import backend as K
Using TensorFlow backend.
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:516: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint8 = np.dtype([("qint8", np.int8, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:517: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_quint8 = np.dtype([("quint8", np.uint8, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:518: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint16 = np.dtype([("qint16", np.int16, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:519: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_quint16 = np.dtype([("quint16", np.uint16, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:520: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint32 = np.dtype([("qint32", np.int32, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:525: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  np_resource = np.dtype([("resource", np.ubyte, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:541: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint8 = np.dtype([("qint8", np.int8, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:542: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_quint8 = np.dtype([("quint8", np.uint8, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:543: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint16 = np.dtype([("qint16", np.int16, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:544: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_quint16 = np.dtype([("quint16", np.uint16, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:545: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  _np_qint32 = np.dtype([("qint32", np.int32, 1)])
/Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:550: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  np_resource = np.dtype([("resource", np.ubyte, 1)])
In [2]:
np.linspace(1, 10, 10)
Out[2]:
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.])
In [3]:
## Some helper functions
def map2layer(x, layer):
    feed_dict = dict(zip([model.layers[0].input], [x.copy()]))
    return K.get_session().run(model.layers[layer].input, feed_dict)

def matplotlib_to_plotly(cmap, pl_entries):
    h = 1.0/(pl_entries-1)
    pl_colorscale = []

    for k in range(pl_entries):
        C = list(map(np.uint8, np.array(cmap(k*h)[:3])*255))
        pl_colorscale.append([k*h, 'rgb'+str((C[0], C[1], C[2]))])
    return pl_colorscale


# Create random data with numpy

def plotly_LSTM_with_outliers(data, error_index_list, time_steps=20):
    fig = go.Figure()
    x_axis = np.linspace(1, time_steps, time_steps)
    for i in range(len(data)):
        data_i = data[i]
        if i in error_index_list:
            fig.add_trace(go.Scatter(x=x_axis, y=data_i,
                        mode='lines',
                        name=i, line=dict(color='red', width=1)))
        else:
            fig.add_trace(go.Scatter(x=x_axis, y=data_i,
                        mode='lines',
                        name=i, line=dict(color='black', width=1)))

    fig.show()
    
    
## some colors
coolwarm_cmap = matplotlib.cm.get_cmap('coolwarm')
coolwarm = matplotlib_to_plotly(coolwarm_cmap, 255)

magma_cmap = matplotlib.cm.get_cmap('magma')
magma = matplotlib_to_plotly(magma_cmap, 255)
In [4]:
# HPCCv1 data
variables_name = pd.read_csv("../data/hpcc_v1/variables_name.csv", header=None)
features = variables_name.values[:,1]


with open("../data/hpcc_v1/X_train_HPCC_1_20.json") as of:
    X_train = np.array(json.load(of))
with open("../data/hpcc_v1/y_train_HPCC_1_20.json") as of:
    y_train = np.array(json.load(of))
with open("../data/hpcc_v1/X_test_HPCC_1_20.json") as of:
    X_test = np.array(json.load(of))
with open("../data/hpcc_v1/y_test_HPCC_1_20.json") as of:
    y_test = np.array(json.load(of))    
In [5]:
## Sort data by target
train_sorted_index = np.argsort(y_train)
X_train = X_train[train_sorted_index]
y_train = sorted(y_train)
In [6]:
# Recreate the exact same model, including its weights and the optimizer
model = tf.keras.models.load_model('../models/HPCCv1_model.h5')

# Show the model architecture
model.summary()
Model: "sequential_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
lstm_1 (LSTM)                (None, 20, 8)             608       
_________________________________________________________________
lstm_2 (LSTM)                (None, 20, 8)             544       
_________________________________________________________________
flatten_1 (Flatten)          (None, 160)               0         
_________________________________________________________________
dense_1 (Dense)              (None, 8)                 1288      
_________________________________________________________________
dense_2 (Dense)              (None, 4)                 36        
_________________________________________________________________
dense_3 (Dense)              (None, 1)                 5         
=================================================================
Total params: 2,481
Trainable params: 2,481
Non-trainable params: 0
_________________________________________________________________
In [7]:
explainer = shap.DeepExplainer(model, X_train)
shap_values_train = explainer.shap_values(X_train)
WARNING:tensorflow:From /Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/shap/explainers/deep/deep_tf.py:502: add_dispatch_support.<locals>.wrapper (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.
Instructions for updating:
Use tf.where in 2.0, which has the same broadcast rule as np.where

SHAP for ALL Instances

In [8]:
################# Plot AVERAGE shap values for ALL observations  #####################
## Consider ABSOLUTE of SHAP values ##
shap.initjs()
shap_average_abs_value_train = np.abs(shap_values_train[0]).mean(axis=0)

x_average_value_train = pd.DataFrame(data=X_train.mean(axis=0), columns = features)
shap.force_plot(0, shap_average_abs_value_train, x_average_value_train)
Out[8]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
In [9]:
shap_values_train_2D = shap_values_train[0].reshape(-1,10)
X_train_2D = X_train.reshape(-1,10)

shap.summary_plot(shap_values_train_2D, X_train_2D, features)

Explain for each instance (input layer)

In [10]:
import seaborn as sns
import matplotlib.pyplot as plt

vmin, vmax = shap_values_train[0].min(), shap_values_train[0].max()

for i, feature in enumerate(features):
    print(feature)

    plt.figure(figsize = (8,6)) 
    tmp = shap_values_train[0][:,:,i].reshape((-1,20))
    print(tmp.shape)
    plot_shap = sns.heatmap(tmp, cmap="coolwarm", vmin= vmin, vmax=vmax)
    plt.show(plot_shap)
    print("-----------")
CPU1 Temp
(312, 20)
-----------
CPU2 Temp
(312, 20)
-----------
Inlet Temp
(312, 20)
-----------
CPU Load
(312, 20)
-----------
Memory Usage
(312, 20)
-----------
Fan1 Speed
(312, 20)
-----------
Fan2 Speed
(312, 20)
-----------
Fan3 Speed
(312, 20)
-----------
Fan4 Speed
(312, 20)
-----------
Power consumption
(312, 20)
-----------

Some Comments

  • CPU1 Temp: Light color at early time steps. It starts bolder from 10th to 18th steps => These steps play an important role in prediction. (recall that output is the sum of 20*10 importance scores)
  • Other features seem contribute slightly on the output.
In [11]:
## Print layers'shape
from keras import backend as K

inp = model.input                                           # input placeholder
outputs = [layer.output for layer in model.layers]          # all layer outputs
functors = [K.function([inp], [out]) for out in outputs]    # evaluation functions

# Testing
layer_outs = [func([X_train]) for func in functors]

for i in range(len(layer_outs)):
    print(layer_outs[i][0].shape)
WARNING:tensorflow:From /Users/chau/opt/anaconda3/envs/py37_env/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:422: The name tf.global_variables is deprecated. Please use tf.compat.v1.global_variables instead.

(312, 20, 8)
(312, 20, 8)
(312, 160)
(312, 8)
(312, 4)
(312, 1)

Explain for LSTM layers

In [12]:
# explain how the input to the ith layer of the model explains the prediction

for i in range(1, 3):
    layer = i
    e = shap.DeepExplainer((model.layers[layer].input, model.layers[-1].output), map2layer(X_train.copy(), layer), )
    shap_values_i, indexes = e.shap_values(map2layer(X_train, layer), ranked_outputs = 1)
    
    vmin, vmax = shap_values_i[0].min(), shap_values_i[0].max()
    print("====================")
    print("===== Layer {} =====".format(i))
    print("shap values shape =", shap_values_i[0].shape)
    for j in range(shap_values_i[0].shape[-1]):
        print("Node {}".format(j+1))
        plt.figure(figsize = (8,6)) 
        tmp = shap_values_i[0][:,:,j].reshape((-1,20))
        print(tmp.shape)
        plot_shap = sns.heatmap(tmp, cmap="coolwarm", vmin= vmin, vmax=vmax)        
        plt.show(plot_shap)
        print("-----------")
    print("-----------------------\n\n") 
====================
===== Layer 1 =====
shap values shape = (312, 20, 8)
Node 1
(312, 20)
-----------
Node 2
(312, 20)
-----------
Node 3
(312, 20)
-----------
Node 4
(312, 20)
-----------
Node 5
(312, 20)
-----------
Node 6
(312, 20)
-----------
Node 7
(312, 20)
-----------
Node 8
(312, 20)
-----------
-----------------------


====================
===== Layer 2 =====
shap values shape = (312, 20, 8)
Node 1
(312, 20)
-----------
Node 2
(312, 20)
-----------
Node 3
(312, 20)
-----------
Node 4
(312, 20)
-----------
Node 5
(312, 20)
-----------
Node 6
(312, 20)
-----------
Node 7
(312, 20)
-----------
Node 8
(312, 20)
-----------
-----------------------


Explain for Dense Layers

In [13]:
for i in range(4,len(layer_outs)):
    layer = i
    e = shap.DeepExplainer((model.layers[layer].input, model.layers[-1].output), map2layer(X_train.copy(), layer), )
    shap_values_i, indexes = e.shap_values(map2layer(X_train, layer), ranked_outputs = 1)
    shap_values_i[0].shape

    print("====================")
    print("====== Layer {} =====".format(i))
    print("shap values shape =", shap_values_i[0].shape)
    for j in range(shap_values_i[0].shape[-1]):
        print("Node {}".format(j+1))
        
        plt.figure(figsize = (8,6)) 
        tmp = shap_values_i[0][:,j].reshape(len(X_train))
        tmp_df = pd.DataFrame(data=[tmp, y_train])
        tmp_df = tmp_df.transpose()
        tmp_df.columns = ["shap_value", "label"]
#         tmp_df.sort_values("label", inplace=True)
        tmp_df["No."] = range(1, len(tmp_df)+1)

        tmp_1 = tmp_df[["shap_value", "No."]].copy()
        tmp_2 = tmp_df[["label", "No."]].copy()

        tmp_1.columns = ["value", "No."]
        tmp_1["type"] = "shape_value"
        tmp_2.columns = ["value", "No."]
        tmp_2["type"] = "label"

        tmp_full = pd.concat([tmp_1, tmp_2])
        
        plot_shap = sns.scatterplot(x="No.", y="value", data=tmp_full, hue="type")
        plt.show(plot_shap)
        print("-----------------------\n\n") 
====================
====== Layer 4 =====
shap values shape = (312, 8)
Node 1
-----------------------


Node 2
-----------------------


Node 3
-----------------------


Node 4
-----------------------


Node 5
-----------------------


Node 6
-----------------------


Node 7
-----------------------


Node 8
-----------------------


====================
====== Layer 5 =====
shap values shape = (312, 4)
Node 1
-----------------------


Node 2
-----------------------


Node 3
-----------------------


Node 4
-----------------------


Prediction vs Label on Train set

In [14]:
y_pred = model.predict(X_train)

tmp_df = pd.DataFrame(data=[y_pred[:,0], y_train])
tmp_df = tmp_df.transpose()
tmp_df.columns = ["prediction", "label"]
tmp_df["No."] = range(1, len(tmp_df)+1)

tmp_df["error (%)"] = np.round (np.abs((tmp_df.label - tmp_df.prediction)*100/tmp_df.label), 2)
tmp_df
Out[14]:
prediction label No. error (%)
0 24.821169 21.0 1 18.20
1 30.603172 30.0 2 2.01
2 32.559383 32.0 3 1.75
3 31.147558 32.0 4 2.66
4 32.798008 33.0 5 0.61
... ... ... ... ...
307 84.002731 85.0 308 1.17
308 84.732216 85.0 309 0.32
309 84.757492 85.0 310 0.29
310 84.751793 85.0 311 0.29
311 85.649818 86.0 312 0.41

312 rows × 4 columns

In [15]:
tmp_df["error (%)"].describe([i*.05 for i in range(20)] + [.97,.98,.985,.99])
Out[15]:
count    312.000000
mean       1.212821
std        1.416160
min        0.000000
0%         0.000000
5%         0.095500
10%        0.211000
15%        0.296500
20%        0.394000
25%        0.510000
30%        0.593000
35%        0.678500
40%        0.760000
45%        0.830000
50%        0.930000
55%        1.030500
60%        1.110000
65%        1.211500
70%        1.300000
75%        1.500000
80%        1.730000
85%        1.917000
90%        2.349000
95%        3.183500
97%        3.366800
98%        4.817800
98.5%      5.490150
99%        5.825900
max       18.200000
Name: error (%), dtype: float64
In [16]:
tmp_1 = tmp_df[["prediction", "No.", "error (%)"]].copy()
tmp_2 = tmp_df[["label", "No.", "error (%)"]].copy()

tmp_1.columns = ["value", "No.", "error (%)"]
tmp_1["type"] = "prediction"
tmp_2.columns = ["value", "No.", "error (%)"]
tmp_2["type"] = "label"

tmp_full = pd.concat([tmp_1, tmp_2])
In [17]:
error_thres = 3
In [18]:
fig = px.scatter(tmp_full, x="No.", y="value", color="type",
                 size='error (%)', hover_data=['error (%)'], opacity=0.6,  
                 color_discrete_map = {"prediction": "blue", "label": "black"})
fig.show()
In [19]:
tmp_full["type"]  = tmp_full.apply(lambda x: "big error" if ((x['error (%)'] > error_thres) & (x['type'] == 'prediction')) else x["type"], 1)
fig = px.scatter(tmp_full, x="No.", y="value", color="type",
                 hover_data=['error (%)'], opacity=0.4,
                 color_discrete_map = {"prediction": "blue", "label": "black", "big error": "red"})
fig.show()
In [20]:
tmp_full.to_csv("./prediction.csv")
tmp_full
Out[20]:
value No. error (%) type
0 24.821169 1 18.20 big error
1 30.603172 2 2.01 prediction
2 32.559383 3 1.75 prediction
3 31.147558 4 2.66 prediction
4 32.798008 5 0.61 prediction
... ... ... ... ...
307 85.000000 308 1.17 label
308 85.000000 309 0.32 label
309 85.000000 310 0.29 label
310 85.000000 311 0.29 label
311 86.000000 312 0.41 label

624 rows × 4 columns

In [21]:
# tmp_full  = pd.read_csv("./prediction.csv")
In [22]:
error_index_list= list(tmp_full[tmp_full.type == "big error"].index)

Show outliers

Input

In [23]:
for i, feature in enumerate(features):
    print(feature)

    plt.figure(figsize = (8,6)) 
    tmp = shap_values_train[0][:,:,i].reshape((-1,20))
    print(tmp.shape)
    plotly_LSTM_with_outliers(tmp, error_index_list)
    print("-----------")
CPU1 Temp
(312, 20)
-----------
CPU2 Temp
(312, 20)
-----------
Inlet Temp
(312, 20)
-----------
CPU Load
(312, 20)
-----------
Memory Usage
(312, 20)
-----------
Fan1 Speed
(312, 20)
-----------
Fan2 Speed
(312, 20)
-----------
Fan3 Speed
(312, 20)
-----------
Fan4 Speed
(312, 20)
-----------
Power consumption
(312, 20)
-----------
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>

LSTM

In [24]:
for i in range(1, 3):
    layer = i
    e = shap.DeepExplainer((model.layers[layer].input, model.layers[-1].output), map2layer(X_train.copy(), layer), )
    shap_values_i, indexes = e.shap_values(map2layer(X_train, layer), ranked_outputs = 1)
    
    print("====================")
    print("===== Layer {} =====".format(i))
    print("shap values shape =", shap_values_i[0].shape)
    for j in range(shap_values_i[0].shape[-1]):
        print("Node {}".format(j+1))
        plt.figure(figsize = (8,6)) 
        tmp = shap_values_i[0][:,:,j].reshape((-1,20))
        print(tmp.shape)
        plotly_LSTM_with_outliers(tmp, error_index_list)
        print("-----------")
    print("-----------------------\n\n") 
====================
===== Layer 1 =====
shap values shape = (312, 20, 8)
Node 1
(312, 20)
-----------
Node 2
(312, 20)
-----------
Node 3
(312, 20)
-----------
Node 4
(312, 20)
-----------
Node 5
(312, 20)
-----------
Node 6
(312, 20)
-----------
Node 7
(312, 20)
-----------
Node 8
(312, 20)
-----------
-----------------------


====================
===== Layer 2 =====
shap values shape = (312, 20, 8)
Node 1
(312, 20)
-----------
Node 2
(312, 20)
-----------
Node 3
(312, 20)
-----------
Node 4
(312, 20)
-----------
Node 5
(312, 20)
-----------
Node 6
(312, 20)
-----------
Node 7
(312, 20)
-----------
Node 8
(312, 20)
-----------
-----------------------


<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>
<Figure size 576x432 with 0 Axes>

Dense

In [25]:
for i in range(4,len(layer_outs)):
    layer = i
    e = shap.DeepExplainer((model.layers[layer].input, model.layers[-1].output), map2layer(X_train.copy(), layer), )
    shap_values_i, indexes = e.shap_values(map2layer(X_train, layer), ranked_outputs = 1)
    shap_values_i[0].shape

    print("====================")
    print("====== Layer {} =====".format(i))
    print("shap values shape =", shap_values_i[0].shape)
    for j in range(shap_values_i[0].shape[-1]):
        print("Node {}".format(j+1))
        
        
        tmp = shap_values_i[0][:,j].reshape(len(X_train))
        tmp_df = pd.DataFrame(data=[tmp, y_train])
        tmp_df = tmp_df.transpose()
        tmp_df.columns = ["shap_value", "label"]
#         tmp_df.sort_values("label", inplace=True)
        tmp_df["No."] = range(1, len(tmp_df)+1)


        tmp_1 = tmp_df[["shap_value", "No."]].copy()
        tmp_2 = tmp_df[["label", "No."]].copy()

        tmp_1.columns = ["value", "No."]
        tmp_1["type"] = "shap_value"
        tmp_2.columns = ["value", "No."]
        tmp_2["type"] = "label"

        tmp_1.loc[error_index_list, "type"]  = "big error"
        tmp_full = pd.concat([tmp_1, tmp_2])

        fig = px.scatter(tmp_full, x="No.", y="value", color="type",
                  opacity=0.6,
                 color_discrete_map = {"shap_value": "blue", "label": "black", "big error": "red"})
        fig.show()

        print("-----------------------\n\n") 
====================
====== Layer 4 =====
shap values shape = (312, 8)
Node 1
-----------------------


Node 2
-----------------------


Node 3
-----------------------


Node 4
-----------------------


Node 5
-----------------------


Node 6
-----------------------


Node 7
-----------------------


Node 8
-----------------------


====================
====== Layer 5 =====
shap values shape = (312, 4)
Node 1
-----------------------


Node 2
-----------------------


Node 3
-----------------------


Node 4
-----------------------


In [26]:
print("Last updated: ", datetime.now())
Last updated:  2020-01-28 10:47:00.965890